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ABSTRACT 


This work studies the use of a wave absorbing control law for vibration 
suppression of flexible spacecraft structures. A major advantage of this method is that it 
does not involve truncation into a finite dimensional mathematical model. A closed loop 
scattering matrix was derived which gives the relationship between incoming waves, 
outgoing waves, sensor, and actuator. The control law was determined by minimizing the 
H-infinity norm of this matrix. The control law was applied to the Naval Postgraduate 
School's Flexible Spacecraft Simulator (FSS) for vibration suppression. The simulator's 
flexible beam was controlled using piezoelectric ceramic wafers as sensors and actuators. 
The H-infinity wave absorbing controller contributed significant damping to the structure, 
especially at the first mode of 1 Hz. Therefore, wave absorbing control and piezoceramic 


sensors and actuators offer a viable approach for vibration suppression of space structures. 
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I. INTRODUCTION 


A. BACKGROUND 


Recent trends in spacecraft design have led to larger and more flexible structures. 
Large structures such as Space Station Freedom pose interesting new problems in 
structural dynamics and control. In addition, modern sensor equipment has become more 
and more accurate, making it highly sensitive to structural vibration. For example, an 
antenna or optical sensor, requiring a high degree of pointing accuracy may be mounted 
ona flexible structure connected to the spacecraft bus. The antenna will be subjected to 
vibration from the spacecraft bus due to slew maneuvers, momentum wheel vibration, and 
thruster firing. Passive vibration isolation can absorb some of the energy. However, 
modern payloads have led to a requirement for active vibration control and isolation. 

There have been two major approaches to vibration control of flexible structures. 
The first, and more conventional approach, is based on modal control of the structure such 
as described by Newman [Ref. 1]. This approach usually suffers from modeling errors and 
state estimation. The use of modal model based control laws has been an area of intense 
research during the last decade. The second approach describes the structural response in 
terms of an elastic disturbance which travels through the structure.[Ref. 2] In this 
approach, compensators are designed to reduce the effects of incoming waves on outgoing 


waves. The advantages of the wave absorbing approach are: 


®Relatively simple to implement 
¢Broadband control 


Does not require finite element analysis or modal model 


B. FOCUS OF THESIS 

This thesis concentrates on three areas. First, a wave absorbing controller will be 
derived using an H-infinity approach similar to one used by Matsuda and Fuji [Ref. 3]. 
Second, the control law will be simulated using a finite element model of the beam to be 
studied. The transfer function of the open loop system will be found using modal 
truncation, then the loop will be closed with the H infinity controller. Last, the control 
law will be implemented on the Naval Postgraduate School's Flexible Spacecraft Simulator 
(FSS). The control law will be applied by a real time control system using Matrix, ™, 
System Build™, and AC-100 by Integrated Systems Incorporated (ISI). The results of 


the experiment will be compared the analytical computer simulation. 


This thesis approaches active control of the Flexible Spacecraft Simulator (FSS) 
using a wave absorbing control theory. The response of the beam 1s viewed as a travelling 
elastic disturbance due to a locally applied force. This approach has been developed by a 
number of sources such as von Flotow and Shafer [Ref. 4], MacMartin and Hall [Ref. 5], 
Matsuda and Fuji [Ref. 3] and Fuji et al. [Ref. 6]. 

1. Dynamic Model 

The FSS setup is shown in Figure 1. The base of the flexible beam 1s fixed to the 
center body with the other end of the beam free. Using the Euler-Bernoulli beam theory 


the dynamics of a flexible beam are given by the following partial differential equation: 


EI + pA = 0 (2.1) 


sensor/actuator 





Il. THEORETICAL ANALYSIS 
A. WAVE ABSORBING CONTROL 
Figure 1. FSS Beam Model 


The boundary conditions are : 


Clamped end: w(0,f)=0, w(0,t) = 0 (2.2) 





Actuator moment: EI — (O,.)=M,. (2.3) 


where w/(x,t) is the lateral displacement, M_, is the bending moment, £/ is the bending 
rigidity, and pA is the mass per unit length. The moment and shear force at the free 
end are zero, but will not be necessary for the derivation. As can be seen in equation 


(2.3), the control forces can be accounted for via the boundary conditions, which 
allows the force term on the right hand side of equation (2.1) to be taken as zero. 
a. Laplace Transform Solution to Beam Equation 


The partial differential equation (2.1) can be analyzed using the 


Laplace transform. The initial conditions, w(x, 0) and w(x, 0) are assumed to be zero. 
The nonzero values of w(x, 0) and w(x, 0) obviously do not make a difference in the 


final result. 


LIwx, )] = Wo, s) = \, w(x, He" dt (2.4) 

L{=*] = s? Wo, 9) (2.5) 
aw, _ [© ow sty, — dt 

Ll = Jo are dt = aa W(x, 5) (2.6) 


The fourth order partial derivative becomes an ordinary derivative since time is 


"frozen" in the s domain. This reduces the Euler-Bernoulli beam equation to: 
d Wr, 
gu + 52W(x, 5) =0 (2.7) 
7 | 


where 





) cee (2.8) 
b. State Space Representation 
The analysis of equation (2.7) can be simplified by using state 


space methods. A state vector containing the cross sectional properties of the beam is 





introduced. 
yi) id 
0 
2 
yal rrtat, (2.9) 
V3 EI 
a 
Y4 WV 
where, 
w= Sw is the lateral velocity 
0 = st is the angular velocity 
M= Ee is the internal bending moment 
see 
O= EI , ~ is the internal shear moment 
x 


The state vector (2.9) can be used to transform equation (2.7) to a set of 


first order state equations in terms of x in the form: 


<= AY (2.10) 


Let us introduce p = = and proceed as follows. 


dy | 


dx =#=§=y) (2.11) 


dy> do dw Md s 
dx ae = ae ~ Rp a¥3 = Py3 (2.12) 


(2.13) 


dy3 _ a dM _ ay _ 
de = Eide = EIS = 4 
Osa _ a £(EI2*) = a 2.14 
ae Ede 7 ae (2.14) 
Using the dynamic system equation (2.7) we can detemine os as 
dtw a: 
a ane (2.15) 
Therefore 
— = —psw = —pw = —py| (2.16) 
This leads to the state equation: 
0 100 
dY 0 0p0 
— = r=AY 2.17 
«~) 9001 _e 
—p000 


The solution of the ordinary differential equation shown above ts given 


by 
¥(x) = e**¥(0) 
= Ue*U'Y(0) (2.18) 


where A is found by similarity transformation of A. This can be done by solving the 


following eigenvalue problem. 


UA =AU 





(2.20) 


I 

hy 

U=| (2.21) 
P 
ry 
Pp 


Slo sleog 
an) ee 
S| slea 


where the eigenvalues are 


l+J 


oy aie (2.22) 


—]+1 
-l-ij 
As can be seen from equations (2.18-22), the solution yields imaginary 
solutions. A more useful solution can be transforming the matrix A or its diagonal 
transformation to real Jordan canonical form. [Ref. 3] As discussed in the next section 
this will also yield a state vector in terms of travelling waves instead of cross sectional 
variables. 
€. State Vector in Terms of Travelling Waves 
The state vector Y from the previous section is defined tn terms of 
the cross sectional variables of the beam. As stated at the beginning of this chapter, we 


desire to study the disturbances in the beam in terms of travelling waves. The matrix A 


can be transformed to real Jordan canonical form by the matrix K . [Ref. 7, pp 144,145] 


Oo © 


CO OO vl wl— 
— mim ot 


The new matrix A can then be formed by the similarity transformation 
A =(K7!U7')A(UK) = P“!AP 


where 


Spe) PSR PP 
7 g Y2p 0 -V2p 
_pi pi pi pi 


Using this transformati-n the state vector can be changed to 


Y=PV 


where the new state vector }’ contains the amplitudes of the travelling wave modes. 


(2.23) 


(2.24) 


(2.25) 


(2.26) 


(2.27) 


This notation is similar to that used in liturature used in microwave circuits. [Ref. 4, pg. 


674] At any point along the beam the motion can be described by four modes, two 


incoming and two outgoing waves as shown in Figure 2. 








Figure 2. Travelling Wave Modes 


The new state vector varies with postion along the beam according to the 


new state equation 


110 0 
dl a Pp 1 | O QO 
—=AV= J5 2.28 
a 7100-1 1 oe 
00-1 -l 
Equation (2.28) has a solution of the form 
Vix) = e* V(0) (2.29) 
The matrix ¢* can be represented by 
1 elix 0 
i | 2.30) 
0 ele ( 


Next the inverse Laplace transform (with respect to x) 
e* =1-'[(s]-A)"'] (2.31) 
is used to find the exponentials along the diagonal of equation (2.30). Let c = E 


then 


S—C 


gl-Ay=| 4 i | (2.32) 
Cc $ 





ei Ae 
(s—c)"+¢ —C §—-C 


Ce ae = C (2.33) 


By taking the inverse Laplace transform using equation (2.31) the first diagonal 


element of equation (2.30) is found. 


re e*“coscx e“simcx 
eo = (2.34) 
Cx - cx 
—e“*sincx e“Coscx 


The second diagonal element can be found in a similar fashion leading to the solution 


e“coscx e“sincx 0 0 
—e“sin cx e* Cos cx 0 0 
V(x) = oe ae V(O) (2.35) 
0) 0 e“coscx e“sincx 
0 0 ~—e~* sin cx e* COS CX 


To determine the wave amplitudes for the pinned and free end of the beam in Figure 


2. the solution is determined for x = L. This leads to the same solution as determined 


by Matsuda and Fuji. [Ref. 3] 


ay e“*coscx e*sincx 0 0 ay 
ay —e“sincx e*coscx 0 0) a 
(L) = * +(0) 


b, 0 0 e“coscx e “sincx by 
b> 0) 0 —e~“*sincx e~* Cos cx b> 


(2.36) 


d. Boundary Conditions 
The boundary conditions are given by equations (2.2) and (2.3). The 


boundary conditions can be transformed to the new coordinates by equation (2.25). 


w 

1000} 0 0 

Peed te — 
a’ 


which can be rewritten as 
2p 0 W2p 0 la 


aoe C) p> ps —p3 ps ar | _ 
00#0 || 0 f2p 0 -2p |} 4 


OF 


es 0 i/2p 0 | ay -| 0 ais 
: M.. : 


This equation can be solved for the outgoing wave amplitudes [5, , >,] as a function of 


the incoming wave modes [ a, ,a,] and the control moment M.. 


b -10 0 |, 
aoe aoe 2 tee (2.40) 
Pp 


Z. Control Law 
a. System Equation in Four Block Form 
The wave equations derived in Section 1 can be written as a 


standard four block form in robust control theory. The equations are as follows. 
b=Sa+ Bu (2.41) 
y=la+Gu (2.42) 


where the matrix 5, which in this special case is called the scattering matrix contains 


the reflection coefficients. The standard four block form is written as 


b SB a 5 
tiered a 


The block diagram for this system is shown in Figure 3. 


y = angular 
velocity 





Figure 3. System Block Diagram 


The controller inputs are of the form 


u= Ky (2.44) 
7M = KO 


Combining compensator K with the output equation given by (2.42) gives 


u=KTa+KGu (2.45) 
OT 

u=Hla (2.46) 
where 

H=K(I-Gky' (2.47) 
or, in other words, 

K = HU+ HG)" (2.48) 


The closed loop system between a and 4 can be determined by combining equation (2.46) 


with (2.41) 


h=Sa+BHTa=(S+BHT)a (2.49) 
Or 
b= S.1a (2.50) 


where S_, is the closed loop scattering matrix. 
For the system derived in Section 1, the four block form equation can be 


determined from the transformation (2.26) and equation (2.40). 


Ww V2p 0 2p 0 |fa 


6 (1)| p> -p p> |) ar 

a wy, (2.51) 
EI ? 0 /2p 0 il) p by 
AV -_pi pi pi pi b> 


The input equation (2.44) can be determined by using the second and third row equations 


in (2.51) above. 
‘ 3 3 3 3 
8 = sp?a\ + 5P?a) — sp2b) + sp?b> (2.52) 
2M. = /2 par - ¥2 pb> (2.53) 


The above two equations are combined producing 


tule 


@ = pia) + p2ar- (2) Mc ali 


Equations (2.40) and (2.54) form an equation appropriate for a control system design. 


h —] Q 0 4 

l oj | 

be eat ee -V2p , 4 a (2.55) 
0 p? Pp? =(4 ” M. 


The first row of this equation can be dropped since 6, cannot be physically controlled by 
M__ This wave mode is called a “far field" mode since it represents a mode transmitted 


along the beam's axis for a long distance. Therefore, the closed loop system Scan be 


determined. 
sa=[0 1 ]+(-v20") Hf p? p? | 
=[-x¥1-Xx | (2.56) 
where 
x= J2Hp? (2.57) 


b. H infinity Compensator Design 

In order to minimize the effect of the incoming waves on the outgoing 
waves, it is desirable to choose a compensator K for which the "magnitude" of the closed 
loop scattering matrix is minimized. This is motivated by drawing energy out of the 
system at the junction point where the actuator is located. 

The magnitude of a matrix can be measured by its norm. One type of 


matrix norm is the H infinity norm which 1s given by 


[IS cill,, = sup(6[S.(7@)]) (2.58) 


© is the largest singular value. The singular values can be found by 
6(S.1) = JAS Sct) (2.59) 


where S., is the complex conjugate of S,;, and A; are the eigenvalues of the matrix in 


parentheses. The eigenvalues are found by 


S'S = oe I ~X(jo) 1—-X(jo) | (2.60) 
K-XGO)X7B) —— X-/OL1-XV)]— | _ gery 
X(j@)[1 — X-yo)] 4- [1 -AVo)][1 - AC) 
Solving the determinant above leads to 
6? =.= 2[ Xa) -$ | Xa) - 4] +4 (2.62) 


The smallest solution for equation (2.62) occurs at X(5) = ; which gives 6 = Bre Now 


equation (2.57) can be solved for H. 


_ 1 yay; 
H= FS) (2.63) 


Using equation (2.48) K can be determined as 
2 a i 
K= 3p 6s)” (2.64) 


As seen in equation (2.44) K gives the relationship between the angular 


velocity, , and the actuator moment, 7. However, the sensor for the FSS beam gives 
the angle 9 for an output. Thus it is desired to find a transfer function from 0 toM. This 


transfer function can be found by taking the Laplace transform of (2.44). 
—M As) = Ks@(s) (2.65) 
Combination of this equation with (2.8) and (2.64) yields 


C(s) = = = Zao : pay Js (2.66) 

This transfer function was derived for a beam clamped at the base and free 

at the end. The transfer function above was also derived for a pinned-free beam by 

Matsuda arid Fuji [Ref. 3] using the H-infinity approach. A similar tranfer function was 

found by MacMartin and Hall [Ref. 5] by minimizing the H infinity norm of the power 

flow into the system. MacMartin and Hal: wiso point out that this solution extracts half of 

the power at all frequencies. This is in contrast to other controllers which are more 

effective at a given mode, but cannot add significant damping at all of the modes at the 
same time. 

The control law given in (2.66) is effectively a "half differentiator". It is 


similar to velocity feedback with a forty five instead of a ninety degree phase lead. In 
order to implement this control law the function /s must be estimated by a rational 





transfer function. As shown by MacMartin and Hall [Ref. 5], this irrational function can 
be approximated over a wide frequency range by a finite number of zeros and poles spaced 
logrithmically along the negative real axis. Figure 4 shows the Bode plot for four zeros 


and four poles. 


(110-4 (10-2 (s410°) (10 


(2.67) 
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Figure 4. Rational Estimate of //S 


B. SENSOR AND ACTUATOR 

The sensors and actuators for the FSS beam are nearly collocated at the base and 
can be seen in Figure 5. Both the sensor and actuator are made of piezoelectric ceramic. 
As seen in the figure, the piezoelectric wafers are stacked. This allows for high sensitivity 


from the sensor and more actuator power in a small area. 
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Figure 5. First Beam Element with Sensor Actuator Pair 


l. Piezoelectric Material 

Piezoelectric ceramics have many advantages as actuators and sensors. Table | 
shows the advantages of using piezoelectric materials in structural vibration control. 
[ Ref. 8] 

The usefulness of piezoelectric ceramics in vibration control is derived from their 


ability to convert electrical energy into mechanical energy and vice versa. When a force 1s 


18 


applied to the material, it creates a voltage proportional to the applied force. Conversly, 


when an electric potential is applied across the ceramic, it creates mechanical strain. 





e  Semor ARO 


high strain sensitivity high stiffness 




















low noise sufficient stress to control vibrations 


low-moderate tempurature sensitivity [good linearity 


easily implemented low power consumption 





easily implemented 






Table 1. Advantages of Piezoelectric Ceramics 


Zi Sensor 
It can be shown that a piezoelectric material produces voltage proportional to 


strain when force 1s applied. [Ref. 9] This voltage is given by the equation 


Ve = ts 3 (1—v)e, (2.68) 


The material constants for the Navy Type I] PZT, which is used in this study, are given in 
Table 2. 







d,, Lateral Strain m/V or Coul/N 1 8e-10 
E Young's Modulus Pascal (N/m’) 6.30e10 
v Poissons Ratio n/a 0.35 
D Abs Permittivity Farad/m or N/V 1.5e-8 






2x 1.905e-4 
1. 588e-3 
4 839e-4 





Sensor Thickness 





Beam Thickness 









Sensor Area 





Table 2. Material Constants for Navy Type IT PZT 


The strain developed at the cross section of the beam at a distance / from the midplane of 


the beam is 


O-w 


Ex = lh (2.69) 
Combine this equation with (2.68) to get 
V(x) =k (2.70) 
where 
Ed 
k = —-ht,—(1 — v) (2.71) 


Since the sensor is distributed along the beam, however, the measured sensor output is 


the average voltage over the length, 1, of the wafer. 


VY, == 0" Ved (2.72) 


V, = +0 (2.73) 


Thus, the voltage output of the sensor is directly proportional to the effective bending 
angle over the sensor element. The value of 7 for the FSS beam is -1.378 x 10° 
rad/volt. 

3. Actuator 

The actuator produces a moment proportional to the voltage across its terminals. 


The electric field across the wafer is related to the voltage by 
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® =— (2.74) 


A lateral strain is induced by this electric field and 1s given by 

€ = da3;D (2.75) 
Furthermore, the stress in the piezoelectric wafer is 

o===Ee, (2.76) 


where A_ is the wafer cross section. If} is the width of the actuator, then equations (2.74) 


to (2.76) can be combined to find the force along the midplane of the actuator. 
Fy, = bEd3V (2.77) 


Since there is a piezoelectric stack on each side of the beam, the moment developed by the 


actuator 1S 





M = 2bEd3\(*)V (2.78) 


Using b = 1.905 cm. and values from Table 2 this equation becomes 
M = (.01675*2)V (2.79) 


Thus, the actuator produces a moment which is proportional to the voltage applied across 


its terminals. 
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Ii. STRUCTURAL RESPONSE SIMULATION 


The open loop transfer function was determined by a modal model which was 
generated by truncating a desired number of modes. The open loop transfer function, G(s), 
and the H-infinity controller, H(s), were then discretized and formed into the closed loop 


system as shown in Figure 6. 


reference pe Beam 
Modal Model 
G(z) 


theta(z) 


H-infinity 
Controller 
H(z) 





Figure 6. System Block Diagram 


A. OPEN LOOP TRANSFER FUNCTION 
The beam dynamics were determined by a modal model. The equation of motion 


for a lumped parameter model is given by the matrix equation 
Mg + Cq+Kq=l() (3.1) 


where 


M is the mass matrix 
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C is the damping matrix 
K is the stiffness matrix 


q is the vector of generalized coordinates 


The beam is modeled with eight fintte elements as shown in Figure 7. The mass 
and stiffness matrices were constructed using the finite element method. [Ref. 10, pp. 
300-328] Element mass and stiffness matrices were calculated and combined to form the 
global mass and stiffness matrices for equation (3.1). These matrices can be seen in the 
program BEAM M in the Appendix A. It should be noted that the vertical dispacements 


and rotation at the base are equal to zero due to clamped boundary conditions. 


element lengths 
e(1) =0.114m lumped mass (0.467 kg) 


e(2-7) = 0.0635 m 


e(8) =0.145m 


HON 


qi! q2 q3 a4 





Figure 7. Finite Element Model 
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The generalized coordinate vector can be decoupled into modal coordinates. [Ref. 11, 


pp. 104-107] 


Wy 


0; 
M=) ° ¢=On=>dn@ (3.2) 


Ws 


Og 


The vector n(t) contains the modal coordinates, and the modal matrix @ contains the 


system eigenvectors normalized with respect to the mass matrix. The number of modes, 
n, for the lumped parameter model 1s equal to the number of generalized coordinates. A 
reasonable solution can be found by truncating the number of modes into a physically 
meaningful number. For example, three modes are often adequate to describe the 
longitudinal vibration of a simple cantilevered beam.[Ref. 11, pp. 292-293] 


Equation (3.2) is substituted for the generalized coordinates in equation (3.1) and 


premultiplied by MD! to get the equation of motion in terms of the modal coordinates. 
Ht Act Aan = FO) (3.3) 


The first coefficient matrix is the modal damping matrix. 


2P1B 00 0 0) 
0 e000 v 
A. =O'CO= 0 Oed OQ (3.4) 
0 O06 0 
0 00 0 2P 16 16 
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The modal stiffness matrix 1s 


ow, 0 00 0 
0a 00 0 
Ay=@O'KDM=| 0 0 € 0 0 (3.5) 
0 0 06¢ QO 
0 000 ow; 
and the modal force vector 1s 
0 
M.- 
0 
F(t)=O'F()=07) * (3.6) 
& 
& 
0 
0 


Note that since the actuator is mounted on the first element, the external moment, M, , is 


applied to the second coordinate, g(2) = 91. 
The modal matix equation (3.3) can be separated into its individual equations in 
the form 
hit 2PiBiNi 5 B; Ni =f u(t) DL Lye (3.7) 
The Laplace transform of this equation ts. 
sy; + 82p;aini(s) + on n:(S) =f; Ui(s) (3.8) 


which leads to the transfer function for n,(s) and U/(s) as follows 
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nits) _ fi 


— 2 
Us) 524295 ;5+0 ; 





(3.9) 


Equation (3.9) gives the transfer function from the input, U/(s), to the 1-th modal 
coordinate,N;, where 7 = 1,2,...,16. 


Now, the system transfer function for a physical coordinate can be determined. 


The output equation for a single input, multiple output system is 


Wt) = C1 q(t) (3.10) 


As discussed in Chapter II, the piezoelectric sensor mounted on the first element measures 


the angle 8;. Thus the output equation becomes 


w)=OH=[010¢0 loa (3.11) 
01(0) = 92(t) (3.12) 


This equation leads to a single input, single output system. Using the transformation to 


modal coordinates given by equation (3.2), the output equation becomes 
i 
81() = & dani) (3.13) 
j= 


where $2, denotes elements of second row of ®. Convert this equation to the s domain 


using the Laplace transform. Therefore, 
H 
0)(s) = > 02 ;Ni(S) (3.14) 
i=| 


Combining this equation with (3.9), and the fact that the input U/(s) -M.(s), yields the open 


loop transfer function. 


pg 


Oi(s) 27 


G(s) = = 3.15 
( ) M-(s) j=] 5242p; a sto? ( ) 


As seen in Figure 8, one way to view equation (3.15) is » single transfer functions 


combined in parallel. This can be easily implemented using the Matlab™ Control Systems 


Toolbox "parallel.m" command. 





Figure 8. Open Loop Transfer Function 


B. DISCRETE TRANSFER FUNCTIONS 
Each transfer function was tranformed to the z (discrete) domain before forming 


the closed loop system. The relationship between the s (continuous) and z (discrete), 1s 
=e! (3.16) 


where T is the sampling period and s = -—a +/ol. 
The open loop transfer function, determined by equation (3.15) was converted to a 


discrete transfer function using a Tustin transformation. The Tustin transformation shown 


below is the Pade approximation to the exponential (3.16). [Ref. 12, pp. 253-282] 
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(3.17) 


The Tustin transformation maintains the frequency response of the continuous system 
while preserving the mapping of the s-plane into the z-plane. 

The transfer function for the estimated H-infinity controller was transformed to the 
z domain by the matched pole zero technique. Recall from Chapter II that the H-infinity 
controller was irrational and was estimated using zeros and poles spaced logrithmically 


along the negative real axis. These zeros and poles can be mapped directly to the z 


domain by replacing each zero or pole term § + @ with its discrete equivalent, Z — ex. 
The gain for the discrete transfer function is chosen so that G(z)le-1 = G(s) 5-0. 
[Ref. 12, pp.304-305] Using Matlab™ the equivalent form of equation (2.64) is 
z+1.000000)(z+.99995 )(z+.995012)(z+.60653 | 
G(z) = jie BOUUULY scam eed Coe) Cae), (3.18) 


(2+.999995)(z+.999500)(z+.95 1229)(z+.006738) 


C. CLOSED LOOP SYSTEM 
Control of the beam was simulated by closing the loop using H infinity controller, 


(,/s ) and derivative (s). Also, a second order filter was incorporated to study the effects 


of filtering the sensor output in the system. The closed loop transfer function was 
determined using the feedback command in Matlab™. 
D. SIMULATION RESULTS 

The open and closed loop Bode plots for the H-infinity controller are shown in 
Figure 9. The plot shows a reduction in magnitude over a broad frequency range as 


expected. 
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Discrete Bode Plot: Open and Closed Loop (Ts=0.01 sec) 
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Figure 9. Open and Closed Loop Bode Plots 
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Figure 10. Forced Response 


The Bode plot becomes inaccurate at frequencies above the Nyquist rate of 314 radians 

per second due to the system being discrete with a sampling period of 0.01 seconds. 
Figure 10 shows the time response to a sinusoidal input of 1 Hz which is the first 

natural frequency of the beam. This plot shows a dramatic reduction in the amplitude of 


the closed loop system for the first two modes (1 and 17 Hz). 
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IV. EXPERIMENTAL ANALYSIS 


A. SET UP 


The wave absorbing controller was implemented on the Flexible Spacecraft 


Simulator at the Naval Postgraduate School. Figure 11 shows an overall picture of the 


system. 





Figure 11. Flexible Spacecraft Simulator 


The spacecraft simulator consists of a center body which floats over a granite table 
on air pads and rotates around an air bearing. A momentum wheel is mounted on the body 
to allow for slew maneuvers and in this experiment is used as a disturbance source. The 


flexible beam is attached to the center body as shown in Figure 12. 
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Figure 12. Flexible Beam 
The beam is controlled by a piezoelectric sensor-actuator pair mounted at the base of the 


beam as seen in Figure 13. The dimensions of the sensor and actuator are shown in Figure 


5, Chapter I. 





Figure 13. Sensor and Actuator 


The digital control system was implemented using Matrix,"”’, System Build’, and 
AC-100'™ control hardware. The control law was designed using block diagrams in 
Or 


System Build™. The system was then compiled, linked and loaded into the AC-100°", a 


real time controller. Figure 14 shows the control software block diagram. 
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Sampling Interval First Sample 
0.0100 0. 


wave mag 
24 ia 
Y= U3*(1 + 01/100*U2) 
2) 
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a] Lo) 


xD 0.018722 +0,03642_+-0.01822 
z*-1.50352 + 0.6564 
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The AC-100™ has an output of plus or minus 10 Volts. This signal was amplified 
by two Kepco™ amplifiers connected in series. This allowed a range of plus or minus 75 
Volts for the piezoelectric actuator input. 

The momentum wheel was used as a disturbance force and is controlled by voltage 
with 1 Volt per 300 RPM. In order to provide a disturbance torque at a desired frequency 
the momentum wheel was spun up to a set RPM, then the wheel RPM was modulated at 
the desired frequency. The control blocks for the momentum wheel are also shown in 
Figure 14. 

B. EXPERIMENTAL RESULTS 

I. Initial Tests 

The initial choice for the digital sampling rate was 100 Hz. Sampling rates of 50, 
200, and 400 Hz. were also used. The 100 and 200 Hz sampling rates appeared to be a 
reasonable choice. 

The main difficulty at the initial stages of the experiment was high frequency noise 
coming from the sensor. Whenever there was a large jump in the sensor signal between 
samples, the output to the actuator would have a large jump in voltage, causing the 
actuator to click. These voltage jumps, possibly as high as 150 Volts, could cause damage 
to the piezoelectic actuator. Thus, different sampling rates were tried. The higher 
sampling rates worked well at reducing the voltage steps between samples, since the 
sampling period was shorter, making the actuator output closer to continous. However, 
the higher sampling rates also allowed more of the higher frequency noise into the system. 
The lower sampling rates, allowed less noise in the system, but had large voltage jumps 
between samples. 

The solution to these problems was to insert a low pass filter at the sensor output. 
A number of filters of different orders and cutoff frequencies were tried. The final choice 
was a second order filter with a damping ratio of 0.707 and a cutoff frequency of 30 
radians per second (4.8 Hz). This filter, combined with a 100 Hz sampling rate, allowed 


good, stable performance for H-infinity and derivative control. 
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2: Free Response of the Beam 

The free response of the beam was tested for the first two modes. The beam was 
set to an initial condition and allowed to vibrate freely. The plots for response of the first 
two modes are shown in Figures 15 and 16. Note that the second mode plot data was 
collected at a rate of 500 Hz for more accuracy. In order to get a good plot of the second 
mode the data was filtered using a high pass FIR filter with a cutoff frequency of 14 Hz. 
This filtered out the dominant first mode and allowed for a better view of the second 
mode. 

The damping ratio of each mode was determined by noting that the amplitude of 


the wave at any time, f, 1s determined by 
A(t) = Ape o®" (4.1) 


where A. is the initial amplitude and @,, is the frequency of that mode. If the amplitude is 
measured at two different points on the plot the damping ratio can be found by solving for 
| A 
In(Zn) 
~ Byl 





(4.2) 


The damping ratios of the first two modes are shown in Table 3. 


Frequency (Hz / (rad/sec) ) Damping Ratio 
0.95 /5.97 


17.0 / 106.8 





Table 3. Free response damping ratios 


The damping ratios determined experimentally were used for the finite element 
program described in Chapter III. The higher modes were estimated to be about the same 
for the program. Also, note that the first two natural frequencies agree well with the 


program frequencies of 1.01 and 17.8 Hz. 
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Figure 15. Free Response, Mode | (0.95 Hz) 
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Figure 16. Free Response, Mode 2 (17 Hz) 
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S Initial Condition Response with Active Control 
In this part of the experiment, the response of the beam to an initial condition was 


tested with two different control laws. The first was the H-infinity (or Js ) 


controller, and the second was a derivative control. As pointed out in Chapter II, the 
H-inifity controller in this case is similar to the derivative control. It provides a phase lead 
of 45 degrees instead of 90 for the derivative controller. Its magnitude rolls off more 
sharply compared to the derivative controller at higher frequencies. 

The results are shown in Figures 17 and 18. A damping ratio for each case was 
calculated in a similar fashion to the free response in the previous section. The measured 


frequencies and calculated damping ratios are shown in Table 4. 


Frequency Damping 
(Hz) Ratio 


0.95 3.9x10° 










Free Response 
H-infinity Control 


Derivative Control 


Table 4. Damping Ratios 


From the results in Table 4 the derivative controller seems to be slightly better than 
the H-infinity controller at this frequency. However, as can be seen in Figure 18, the 
derivative controller began to excite the second mode after about 30 seconds. The 
derivative controller was also much more sensitive to noise in the system and often caused 
spikes in the voltage signal sent to the actuator. The H-infinty controller, on the other 
hand. was less sensitive to noise and gave a smoother input to the actuator. The H-infinity 


controller turns out to produce reliable performance and was more stable. 
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Figure 17. Initial Condition Response using H-infinity Control 





0 500 1000 4500 2000 2500 3000 3500 4000 
N {secx100) 
Beam Angle — Derivative Control (T3.01) 


Figure 18. Initial Condition Response using Derivative Control 
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4. Forced Response 

The forced response of the closed loop system was tested using the momentum 
wheel as a disturbance source. Different runs were made using different frequencies. 
Both H-infinity and derivative control laws were implemented. 

As discussed previously, the system had stability problems with high frequency 
noise and large voltage jumps. Thus, there was a necessity to use a low pass filter. 
Initially first and second order filters with 30, 50 and 100 rad/sec cutoff frequencies were 
used. The most stable performance was achieved using a second order 30 rad/sec filter. 
This reduced the effective bandwidth of the controller to 4.8 Hz. Also, the filter began to 
reduce the amount of phase lead of the controller as the disturbance frequency came near 
4.8 Hz. 

The controller did, however, work well within the low pass bandwidth with 
significant active control action near the first natural mode of 1 Hz. The results of 
closing the control loop are given in Figures 19 and 20. The decrease in amplitude for 


both controllers was calculated and is shown in Table 5. 


(Initial Final) amplitude | Decrease in amplitude(dB) 


H-infinity 


Derivative 





Table 5. Controller performance (1 Hz disturbance) 


The H-infinity controller performed better in this case. A significant reason for 
this was that the H-infinty controller was less sensitive to noise and, therefore, the gain 
could be increased to get full use of the 150 V peak to peak output of the actuator 
amplifier. The derivative controller became unstable if the gain was increased too much, 
even with the low pass filter. Figure 20 shows that the derivative controller begins to 
excite higher modes after about 55 seconds. The H-infinity controller, combined with the 
second order filter, showed no tendency to excite the higher modes when subjected to the 


disturbance. The H-infinity controller had the best performance with a forced response. 
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Figure 19. Forced Response, H-infinity Control (1 Hz disturbance) 
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Figure 20. Forced Response, Derivative Control (1 Hz disturbance) 
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V. SUMMARY AND CONCLUSIONS 


Active control of flexible spacecraft structures has become an important topic. 
Larger and lighter structures supporting more sensitive equipment have led to the need for 
active control of spacecraft structures. This thesis studied a wave absorbing approach to 
this control problem. A closed loop scattering matrix was derived which determines the 
relationship between incoming waves, outgoing waves, sensor and actuator. The control 
law was found by minimizing the magnitude of this closed loop scattering matrix. 
Theoretically, this approach allows for a relatively simple control law, broadband control, 
and does not require a modal model of the system for operation. 

This control law was then applied to the Flexible Spacecraft Simulator (FSS). The 
sensing and actuation for the flexible beam was implemented by piezoelectric ceramic 
waters mounted at the base of the beam. The control law was tested by simulation on a 
computer, and then by applying it to the hardware in the lab. 

The computer simulation showed that the controller was useful over a broad range 
of frequencies below the Nyquist rate of the digital system. The results from the FSS 
confirmed this for frequencies below 5 Hz. At higher frequencies the controller had some 
problems due to sensor noise requiring a low pass filter, and thus limiting the bandwidth of 
the system. At the first natural frequency of | Hz the H-infinity controller performed well 
and increased the damping ratio of the system from 3.93 x 10° to 3.7 x 10°. The 
controller also decreased the amplitude due to a 1 Hz sinusoidal disturbance by 14 dB. 

Some recommendations for further study include adding a weighting function to 
the terms in the closed loop scattering matrix. This would allow for weighting the output 
of the controller at desired frequencies. It also allows for implementation of a low pass 
filter before the H-infinity norm of the matrix is minimized. Another area of study would 
be the use of analog filters on either the sensor output or the actuator input. The former 
could be used to filter out noise from the sensor. The latter would be useful in smoothing 


out the signal sent from the power amplifiers to the piezoelectric actuators. The filters, 
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however, will change the phase of the closed loop system, which is another area to be 
addressed. Lastly, the sensor and actuator were assumed to be collocated for this study. 
Their centers are 1.96 inches apart. This should not be a problem for the lower modes of 
the system. However, it is possible that at higher frequencies, the phase difference 
between the wave at the sensor and actuator could cause a stability problem. This could 
be alleviated by solving the scattering matrix for a noncollocated sensor and actuator. 

In conclusion, the use of a wave absorbing controller with piezoelectric sensors 
and actuators could be a viable approach to vibration suppression of space structures. As 
discussed by von Flowtow [Ref. 13] and Fuji, Ohtsuka and Murayama (Ref. 6], 
eventually, this approach can be applied to large flexible spacecraft structures with 
multiple noncollocated sensors and actuators. The wave absorbing control approach 
combined with new sensor and actuator technology such as piezoelectric ceramics 


presents a wide area for further research. 
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RII O II IGG DIU ISI GICIGIGCI GI IGIOICISIO ISO II TI III TAR IA IR IIIS INI LES 


Beam Simulator * 

* 

This program models a beam 27 inches long with lumped masses* 
to change the natural vibration frequencies. The beam is * 
modeled using a finite element model with 8 elements. The * 
transfer function (theta/moment at element 1) of the beam is * 
found using modal truncation. The program allows a variety of * 
control laws to be tested by the user. Included are the H * 
infinity controller (sqrt(s)), Derivative, and Integral * 
control. Also a second order filter with a variable cutoff * 
frequency can be applied to the H infinity or Derivative * 
* 

* 

* 


CL © ad? oi ad od al? lO alP ol? al? oP ol? ad? al? oP 


controller. 
IG ICICI GIDC DIGI SIGIOICICISIOSIGGICICICICICIO CISCO TTT ITO TI IIE LIES 
lear 
m=.112014; - linear mass density (kg/m) 
ml=.4672; * lumped mass (kg) 
T1=2.0436e-4; ¢ lumped mass rotary moment of inertia (kgm 2) 
h=.0635; ¢ elemental length (m) 
h1i=.1143,; % Ist element length (m) 
h8=.14478; 2 8th element length (m) 
IT=8.46825e-12; % beam cross sectional moment of inetria (kgm 2) 
E=71e9; % Youngs modulus (aluminum) 
n=8; > # of elements 
dampf=0 .0036; ° Estimated damping ratio for beam (Exp from ist 2 
modes ) 
f=zeros(2*n,1); % force vector (actuator) 
£(2)=1; 
kh=0.182926; ° H infinity gain (due to beam properties) 
T=0.01; ¢ Sampling period for discretized system 


° Elemental mass and stiffness 


mill=m*h*[156 22*h;22*h 4*h*2)/420; 
m22=m*h*(156 ~22*h;-22*h 4*h°2)]/420; 
mi2=m*h*(54 -13*h;13*h —3*h°2)/420; 
m21=m12'; 

mileg=m*h8*[156 22*h8;22*h8 4*h8°~2]/420; 
m22e8=m*h8*{156 -22*h8;—-22*h8 4*h8°2]/420; 
m12e8=m*h8*(54 —13*h8;13*h8 —-3*h8°2)]/420; 
m21e8=m12e8'; 

Global mass matrix 


°o 


gm(1:2,1:2)=m22+m11; 
gm(1:2,3:4)=m12; 
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for i=1:n-2 
gm(2*i+1:2*it2,2*i-1:2*1)=m21; 
gm(2*i+1:2*i+2,2*1+1:2*1+2)=m22+m11; 
gm(2*i+1:2*i1+2,2*1+3:2*1+4)=m12; 

end 


gm(2*n-1:2*n,2*n-3:2*n-2)=m21; 
gm(2*n-1:2*n,2*n-1:2*n)=m22; 


gm(13:14,13:14)=m22+m11e8; 
gm(15:16,13:14)=m21e8; 
gm(13:14,15:16)=m12e8; 
gm(15:16,15:16)=m22e8; 


= Add lumped mass (possible 9 nodes) 
MI=[{ml 0; 0 I1]; 

node=8; 

p=2*(node-1); 
gm(p-1:p,p-1:p)=gm(p-1:p,p-1:p)*MI; 


2 Flemental stiffness matrices 


kK11=E*1I*(12 6*h;6*h 4*h*2]/h*3; 
k22=E*1I*(12 -6*h;-6*h 4*h°2)/h° 3; 
k12=E*I*[-12 6*h;-6*h 2*h°2)/h 3; 
k21=k12'; 


k11le8=E*1I*(12 6*h8;6*h8 4*h8°2)/h8°3; 
k22e8=E*I*([12 -6*h8;-6*h8 4*h8°2]/h8°3; 
k12e8=E*I*[(-12 6*h8;-6*h8 2*h8°2)]/h8° 3; 
k21e8=k12e8'; 


% global stiffness matrix 
gk(1:2,1:2)=k22+k11; 
gk(1:2,3:4)=k12; 


for i=l:n-2 
gk (2*it+1:2*i+2,2*1i-1:2*1)=k21; 
gk (2*i+1:2*i+2,2*i+1:2*i+2)=k22+k11; 
gk(2*it1:2*i+2,2*1+3:2*1i+4)=k12; 

end 


gk(2*n-1:2*n,2*n-3:2*n-2)=k21; 
gk(2*n-1:2*n,2*n-1:2*n)=k22; 


gk(13:14,13:14)=k22+k11e8; 
gk(15:16,13:14)=k21e8; 
gk(13:14,15:16)=k12e8; 
gk(15:16,15:16)=k22e8; 
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% Find mode shapes and natural frequencies 


[V,D]=elg(gk, gm); 
[lambda,k]}=sort(diag(D)); 
V=V(:,k); 


factor=diag(V'*gm*V); 
phi=V*inv(sqrt(diag(factor) )); 
omega=diag(sqrt(phi'*gk*phi) ); 
freq=omega/2/pi; 

modf=phi'*f; 


% Calculate the open loop transfer function of the beam at 
% each frequency for the first md modes (truncate). 
md=input('Number of modes (1<#<8): '); 

s2=diag(ones(md) ); 

s1=2*dampf*omega(1:md); 

sO=comega(1:md).*omega(1:md); 


num=phi(2,1:md)'.*modf(1:md); % sensor measures theta # 2 
den=[s2 sl sO]; 
[numg,deng]=parallel(num(1),den(1,1:3),num(2),den(2,1:3)); 
for i=3:md 

(numg, deng]=parallel(numg,deng,num(i),den(i,1:3)); 
end 


% Discretize open loop system 
(numgd, dengd)=c2dm(numg,deng,T, 'tustin'); 


% Control Transfer Function 
flagl='y'; 
while flagl=='y' 


disp('(1) Integral') 

disp('(2) Derivative & Filter') 
disp('(3) sqrt(s)') 

disp('(4) sqrt(s) & Filter: ') 
case=input('Enter Controller Type: '); 
gain=input('Enter gain factor: '); 


Estimate square root of s 
zrs = [-.0001 -.01 -1 -100]; 
pls = [--001 -.1 -10 -1000]; 
nums=poly(zrs); 
dens=poly(pls); 
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(numsd, densd]=c2dm(nums,dens,T, 'matched' ); 


Noise filter 

we=input('Enter Filter cutoff frequency (rad/sec): '); 
dmpf=.707; 

numf=wc’ 2; 

denf=(1 2*dmpf*wc wc'2]; 
(numfd,denfd]=c2dm(numf,denf,T,'tustin'); 


oe 


if case == 

% Integral control 
numhd=gain*[T 0]; 
denhd=[1 -1]; 

elseif case == 
% Derivative 
numhd=conv(numfd, {1 -1]); 
denhd=conv(denfd,[T 0]); 

elseif case == 
numhd=gain*numsd,; 
denhd=densd; 

elseif case >= 4 
(numhd, denhd}]=series(numfd,denfd,gain*numsd,densd) ; 


end 


¢ Closed loop system 
rnumcld, dencla]=feedback (numgd, dengd, numhd,denhd,~-1); 


% Frequency Domain Evalutaion. 
< Plot Bode for open and closed loop systems. 


w=logspace(0,10g10(600),200); 
[mg,pg]=dbode(numgd,dengd,T,w); 
[mc , pc ]=dbode(numcld,dencld,T,w); 


figure(1) 
semilogx(w, 20*log10(mg) ,w,20*10g10(mc) ) 
grid 
title('Discrete Bode Plot: Open and Closed Loop (Ts=0.01 sec)') 
xlabel('omega (rad/sec)') 
ylabel('|G(s){') 


figure(2) 

semilogx(w,pg,W,pc) 

grid 

title('Discrete Bode Plot: Open and Closed Loop (Ts=0.01 sec)') 
xlabel('omega (rad/sec)') 

ylabel('Phase' ) 


Time Domain Evaluation 
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9° 


¢ Plot output of open and closed loop system. 


flag='y'; 
while flag=='y' 


dist=input('Enter disturbance frequency (HZ): '); 
tmax=20/dist; 

nt=0:T:tmax,; 

u=0.05*sin(dist.*6.28*nt) ; 

(yg, xg]=dlsim(numgd,dengd,u); 
(yel,xcl]=dlsim(numcld,denclid,u); 


figure(3) 

plot(nt,yg,nt,ycl) 

grid 

title('Response of System to Sinusoidal Input (Ts=0.01 sec)'); 
Xlabel('time (seconds)'); 

ylabel('theta (rad)'); 


flag=input('Try another disturbance frequency? (y/n) bors! 
end 


flagl=input('Try another control law / filter?(y/n) ','S'); 
end 
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Natural Frequencies 


1.0e+003 * 


0.00101324976094 
0.01775301423174 
0.04638342775662 
0.07923661265561 
0.11935988578432 
0.22849778748830 
0.378 13429140721 
0.56949274733975 
0.61770877977266 
0.79596 117358377 
1.16976439767046 
1.54850209877084 
2.04634252642825 
2.68 166333705337 
3.4539261458164] 
42272623 1062265 


APPENDIX B 


Finite Element Model 


5] 


1.0e+003 * 
Columns | through 7 


0.0000 
0.0013 
0.0002 
0.0023 
0.0003 
0.0032 
0.0006 
0.0039 
0.0008 
0.0044 
0.0011 
0.0046 
0.0014 
0.0048 
0.0021 
0.0048 


0.0005 
0.0142 
0.0016 
0.0195 
0.0028 
0.0165 
0.0036 
0.0066 
0.0035 
-0.0085 
0.0024 
-0.0266 
0.000] 
-0.0458 
-0.0072 
-0.04520 


0.0015 
0.0412 
0.0043 
0.04) 1 
0.0061 
0.0118 
0.0057 
-0.0255 
0.0032 
-(). 0469 
0 0004 
-) 0343 
-0.0002 
00197 
0 0062 
0.0540 


0.0008 
0.0208 
0.0020 
0.0108 
0.0019 
-0.0131 
0.0005 
-0.0279 
-0.0011 
-0.0193 
-0.0015 
00093 
0.0001 
0.0422 
-0.0126 
-() [382 


0.0035 
0.0782 
0.0065 
-0.0003 
0.0030 
-0.1000 
-0.0038 
-0.0919 
-0.0065 
00141 
-0.0028 
0 0800 
0.0001 
-0.0198 
0.0016 
0.0223 


0.0054 
0.0878 
0.0046 
-0.1177 
-0.0045 
-0.1100 
-0.0043 
0.1153 
0.0048 
0.1126 
0.0053 
-0.0923 
-0.0001 
0.0084 
-0.0003 
-0.0063 


0.0067 
0.0479 
-0.0007 
-0.2064 
-0.0054 
0.1106 
0.0055 
0.1064 
0.0005 
-0 2069 
-0 0067 
0.0524 
0.0001 
-0.0049 
0.0001 
0.0054 


0.0070 
-0.0452 
-0.0056 
-0.1344 

0.0022 

0.2379 

0.0021 
-().2392 
-0.0055 

0.1376 

0.0070 

0.0415 
-0.0001 

0.0031 
-0.0004 
-O.OTR6 


0.0001 
-0.0011 
-0.0001 
-0.0015 
0.0001 
0.0033 
0.0000 
-0.0042 
-0.0001 
0.0037 
0.0001 
-0.0013 
0.0001 
0.0079 
0.0221 
1 1674 


Mode Shapes 


-0.0058 
0.1629 
0.0063 

-0.0945 

-0.0067 
0.0333 
0.0067 
0.0306 

-0.0063 

-0.0920 
0.0058 
0.1607 
0.0000 
0.0022 
0.0001 
0 0037 


().0038 
-0.3549 
-0.0023 
0.4203 
0.0008 
-0.4591 
0.0008 
0.4587 
-0.0023 
-(). 4191 
0.0038 
0.3530 
0.0000 
0.0015 
0.0000 
0.0009 


-0.0018 
0.5559 
-0.0028 
-0.5095 
0.0058 
0.2047 
-(1).0058 
0.2072 
0.0028 
-0.5107 
0.0019 
0.5551 
0.0000 
0.0013 
0.0000 
0.0005 


-0.0007 
-0.7169 
0.0060 
0.1693 
-0.0035 
0.5899 
-0.0035 
-0. 5885 
0.0060 
-0.1717 
-0.0006 
0.7171 
0.0000 
0.0011 
0.0000 
0.0003 


0.0027 
0.8024 
-0.0048 
0.5259 
-0 0039 
-0.6777 
0.0039 
-0.679] 
0.0048 
0.5242 
-0.0027 
0.8034 
0.0000 
0.0008 
0.0000 
0.0001 


-().0036 
-0 7484 
0.0005 
-1.1435 
0.0043 
-0.5395 
0.0043 
0.5380 
0.0005 
1 1433 
-0.0036 
0.7497 
0.0000 
0.0006 
0 0000 
0.0001 


0.0027 
0.4795 
0.0019 
1 0360 
0.0007 
1.3528 
-0.0007 
1.3530 
-0.0019 
1.0367 
-0.0027 
0 4805 
0.0000 
0.0003 
0.0000 
0.0000 
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